//
// File: interp1.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 20-Mar-2025 09:03:00
//

// Include Files
#include "interp1.h"
#include "rt_nonfinite.h"
#include "coder_array.h"
#include "omp.h"
#include <cmath>

// Type Definitions
struct struct_T {
  coder::array<double, 2U> breaks;
  coder::array<double, 2U> coefs;
};

// Function Definitions
//
// Arguments    : const ::coder::array<double, 2U> &y
//                const ::coder::array<double, 2U> &xi
//                ::coder::array<double, 2U> &yi
//                const ::coder::array<double, 2U> &x
// Return Type  : void
//
namespace coder {
void interp1SplineOrPCHIP(const ::coder::array<double, 2U> &y,
                          const ::coder::array<double, 2U> &xi,
                          ::coder::array<double, 2U> &yi,
                          const ::coder::array<double, 2U> &x)
{
  array<double, 2U> dvdf;
  array<double, 2U> dx;
  array<double, 2U> md;
  array<double, 2U> s;
  struct_T pp;
  double v;
  double xloc;
  int high_i;
  int low_i;
  int low_ip1;
  int mid_i;
  int numTerms;
  int szs_idx_1;
  boolean_T has_endslopes;
  has_endslopes = (y.size(1) == x.size(1) + 2);
  if ((x.size(1) <= 2) || ((x.size(1) <= 3) && (!has_endslopes))) {
    int yoffset;
    has_endslopes = (y.size(1) == x.size(1) + 2);
    if (x.size(1) <= 2) {
      if (has_endslopes) {
        yoffset = 4;
      } else {
        yoffset = 2;
      }
    } else {
      yoffset = 3;
    }
    pp.coefs.set_size(1, yoffset);
    if (x.size(1) <= 2) {
      if (has_endslopes) {
        double pp1_coefs[4];
        double dnnm2;
        dnnm2 = y[0];
        if (x.size(1) - 2 >= 0) {
          double dzzdx;
          double r;
          double szdvdf_idx_1;
          r = x[1] - x[0];
          szdvdf_idx_1 = (y[2] - y[1]) / r;
          dzzdx = (szdvdf_idx_1 - dnnm2) / r;
          szdvdf_idx_1 = (y[y.size(1) - 1] - szdvdf_idx_1) / r;
          pp1_coefs[0] = (szdvdf_idx_1 - dzzdx) / r;
          pp1_coefs[1] = 2.0 * dzzdx - szdvdf_idx_1;
          pp1_coefs[2] = dnnm2;
          pp1_coefs[3] = y[1];
        }
        szs_idx_1 = pp.coefs.size(1);
        pp.coefs.set_size(1, pp.coefs.size(1));
        for (int nxm1{0}; nxm1 < szs_idx_1; nxm1++) {
          pp.coefs[nxm1] = pp1_coefs[nxm1];
        }
      } else {
        pp.coefs[0] = (y[1] - y[0]) / (x[1] - x[0]);
        pp.coefs[1] = y[0];
      }
      pp.breaks.set_size(1, x.size(1));
      szs_idx_1 = x.size(1);
      for (int nxm1{0}; nxm1 < szs_idx_1; nxm1++) {
        pp.breaks[nxm1] = x[nxm1];
      }
    } else {
      double dnnm2;
      double szdvdf_idx_1;
      szdvdf_idx_1 = x[1] - x[0];
      dnnm2 = (y[1] - y[0]) / szdvdf_idx_1;
      pp.coefs[0] = ((y[2] - y[1]) / (x[2] - x[1]) - dnnm2) / (x[2] - x[0]);
      pp.coefs[1] = dnnm2 - pp.coefs[0] * szdvdf_idx_1;
      pp.coefs[2] = y[0];
      pp.breaks.set_size(1, 2);
      pp.breaks[0] = x[0];
      pp.breaks[1] = x[2];
    }
  } else {
    double dnnm2;
    double r;
    double szdvdf_idx_1;
    int nxm1;
    int yoffset;
    nxm1 = x.size(1) - 1;
    if (has_endslopes) {
      szdvdf_idx_1 = static_cast<double>(y.size(1)) - 3.0;
      szs_idx_1 = y.size(1) - 2;
      yoffset = 1;
    } else {
      szdvdf_idx_1 = static_cast<double>(y.size(1)) - 1.0;
      szs_idx_1 = y.size(1);
      yoffset = 0;
    }
    s.set_size(1, szs_idx_1);
    dvdf.set_size(1, static_cast<int>(szdvdf_idx_1));
    dx.set_size(1, x.size(1) - 1);
    for (int k{0}; k < nxm1; k++) {
      szdvdf_idx_1 = x[k + 1] - x[k];
      dx[k] = szdvdf_idx_1;
      szs_idx_1 = yoffset + k;
      dvdf[k] = (y[szs_idx_1 + 1] - y[szs_idx_1]) / szdvdf_idx_1;
    }
    for (int k{2}; k <= nxm1; k++) {
      s[k - 1] = 3.0 * (dx[k - 1] * dvdf[k - 2] + dx[k - 2] * dvdf[k - 1]);
    }
    if (has_endslopes) {
      szdvdf_idx_1 = 0.0;
      dnnm2 = 0.0;
      s[0] = y[0] * dx[1];
      s[x.size(1) - 1] = dx[x.size(1) - 3] * y[x.size(1) + 1];
    } else {
      szdvdf_idx_1 = x[2] - x[0];
      dnnm2 = x[x.size(1) - 1] - x[x.size(1) - 3];
      s[0] = ((dx[0] + 2.0 * szdvdf_idx_1) * dx[1] * dvdf[0] +
              dx[0] * dx[0] * dvdf[1]) /
             szdvdf_idx_1;
      s[x.size(1) - 1] =
          ((dx[x.size(1) - 2] + 2.0 * dnnm2) * dx[x.size(1) - 3] *
               dvdf[x.size(1) - 2] +
           dx[x.size(1) - 2] * dx[x.size(1) - 2] * dvdf[x.size(1) - 3]) /
          dnnm2;
    }
    md.set_size(1, x.size(1));
    md[0] = dx[1];
    md[x.size(1) - 1] = dx[x.size(1) - 3];
    for (int k{2}; k <= nxm1; k++) {
      md[k - 1] = 2.0 * (dx[k - 1] + dx[k - 2]);
    }
    r = dx[1] / md[0];
    md[1] = md[1] - r * szdvdf_idx_1;
    s[1] = s[1] - r * s[0];
    for (int k{3}; k <= nxm1; k++) {
      r = dx[k - 1] / md[k - 2];
      md[k - 1] = md[k - 1] - r * dx[k - 3];
      s[k - 1] = s[k - 1] - r * s[k - 2];
    }
    r = dnnm2 / md[x.size(1) - 2];
    md[x.size(1) - 1] = md[x.size(1) - 1] - r * dx[x.size(1) - 3];
    s[x.size(1) - 1] = s[x.size(1) - 1] - r * s[x.size(1) - 2];
    s[x.size(1) - 1] = s[x.size(1) - 1] / md[x.size(1) - 1];
    for (int k{nxm1}; k >= 2; k--) {
      s[k - 1] = (s[k - 1] - dx[k - 2] * s[k]) / md[k - 1];
    }
    s[0] = (s[0] - szdvdf_idx_1 * s[1]) / md[0];
    nxm1 = x.size(1);
    szs_idx_1 = s.size(1) - 1;
    pp.coefs.set_size(s.size(1) - 1, 4);
    for (int k{0}; k <= nxm1 - 2; k++) {
      double dzzdx;
      szdvdf_idx_1 = dvdf[k];
      dnnm2 = s[k];
      r = dx[k];
      dzzdx = (szdvdf_idx_1 - dnnm2) / r;
      szdvdf_idx_1 = (s[k + 1] - szdvdf_idx_1) / r;
      pp.coefs[k] = (szdvdf_idx_1 - dzzdx) / r;
      pp.coefs[szs_idx_1 + k] = 2.0 * dzzdx - szdvdf_idx_1;
      pp.coefs[(szs_idx_1 << 1) + k] = dnnm2;
      pp.coefs[3 * szs_idx_1 + k] = y[yoffset + k];
    }
    pp.breaks.set_size(1, x.size(1));
    szs_idx_1 = x.size(1);
    for (nxm1 = 0; nxm1 < szs_idx_1; nxm1++) {
      pp.breaks[nxm1] = x[nxm1];
    }
  }
  szs_idx_1 = xi.size(1) - 1;
#pragma omp parallel for num_threads(omp_get_max_threads()) private(           \
    v, numTerms, high_i, low_i, low_ip1, xloc, mid_i)

  for (int b_k = 0; b_k <= szs_idx_1; b_k++) {
    v = xi[b_k];
    if (std::isnan(v)) {
      yi[b_k] = rtNaN;
    } else {
      numTerms = pp.coefs.size(1);
      if (!std::isnan(v)) {
        high_i = pp.breaks.size(1);
        low_i = 1;
        low_ip1 = 2;
        while (high_i > low_ip1) {
          mid_i = (low_i >> 1) + (high_i >> 1);
          if (((low_i & 1) == 1) && ((high_i & 1) == 1)) {
            mid_i++;
          }
          if (xi[b_k] >= pp.breaks[mid_i - 1]) {
            low_i = mid_i;
            low_ip1 = mid_i + 1;
          } else {
            high_i = mid_i;
          }
        }
        xloc = xi[b_k] - pp.breaks[low_i - 1];
        v = pp.coefs[low_i - 1];
        for (high_i = 2; high_i <= numTerms; high_i++) {
          v = xloc * v +
              pp.coefs[(low_i + (high_i - 1) * (pp.breaks.size(1) - 1)) - 1];
        }
      }
      yi[b_k] = v;
    }
  }
}

} // namespace coder

//
// File trailer for interp1.cpp
//
// [EOF]
//
